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Abstract 

How do single cell fate decisions induced by activation of key signaling proteins above threshold 
• concentrations within a time interval are affected by stochastic fluctuations in biochemical reac- 

o 

tions? We address this question using minimal models of stochastic chemical reactions commonly 
found in cell signaling and gene regulatory systems. Employing exact solutions and semi-analytical 
methods we calculate distributions of the maximum value (N) of activated species concentrations 
(Pmax(N)) and the time (t) taken to reach the maximum value (P ma x(t)) within a time window in 
the minimal models. We find, the presence of positive feedback interactions make P ma x(N) more 
spread out with a higher "peakedness" in P ma x(t). Thus positive feedback interactions may help 
single cells to respond sensitively to a stimulus when cell decision processes require upregulation 
of activated forms of key proteins to a threshold number within a time window. 
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I. INTRODUCTION 



Decisions made at the single cell level enable organisms to respond to changes in the local 
environment. Such decisions are usually processed upon upregulation of specific proteins, 
transcription factors, or soluble molecules that help cells to communicate with each other. 
These activation events often require concentrations of few key proteins to reach a threshold 
level within a time window. Examples of such responses include, activation of immune 
cells such as T cells triggered by a threshold number of pathogenic peptides all or 
none maturation of oocytes in the frog Xenopus laevis induced by different concentrations 
of progesterone [2j, or switch like activation of Lac genes regulating lactose metabolism 
produced by a threshold concentration of stimulus in E. coli [3j . 

However, every cell in a cell population interacting with stimuli possesses unique temporal 
profiles of concentrations of activated signaling molecules or genes. This cell to cell vari- 
ability in the kinetics occurs due to the inherent stochastic nature of associated biochemical 
processes (or intrinsic noise) [4HSI [8] and variations in expression levels of genes and pro- 
teins (or extrinsic noise) [8]. Therefore, the threshold for activation for a specific signaling 
molecule and the time window within which the signaling molecule should be activated to 
influence cell functions can change from cell to cell [2]. How do nonlinearities commonly 
found in biochemical signaling networks, such as positive feedbacks, help cells to respond 
to these variations? We address this question in the article, in particular, we investigate 
the role of positive feedback interactions which are often responsible for producing all or 
none responses in signaling or gene regulatory kinetics. We use a minimal model for a linear 
and a positive feedback interaction in a simple chemical reaction describing activation of a 
single chemical species representing a key signaling protein or a gene. Since many positive 
feedback interactions [U El [7] can be reduced to this form the results from the model will be 
relevant for a wide range of biological systems. 

We consider a minimal biochemical process, C^C*, describing production and deacti- 
vation of the activated species C* which needs to reach a threshold concentration (say N) 
in a time interval [0, T] in order to mediate a functional response. Due to the stochastic 
fluctuations in the kinetics, the threshold concentration of C* could occur at different times 
(Fig. 1) or even stay below the threshold level in the time window. Therefore, knowing the 
distribution of the number (n) of molecules of C* at a time T will not reveal if the concen- 
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tration of C* attained the threshold level at an earlier time. However, knowledge of the joint 
distribution of the maximum number (A) of C* and the time t (0 < t < T) when this value 
was attained in a temporal profile describing the kinetics of C* in a single cell will inform 
us if the cell was able to cross the threshold in the time interval [0,T]. Such distributions 
are regularly dealt with in Extreme Value Theory (EVT) where extreme value distributions 
for identically distributed independent random variables have been studied extensively [9]. 
Analysis of extreme value distributions for correlated random variables has been a topic of 
intense research in the recent years due to its application in physics fTUl [TT1 [T8l fT9l [24] . climate 
science [25] . finance [15]. and, population [T3l fT6] and cell biology [14J. Application of such 
distributions in stochastic biochemical reaction kinetics has been initiated only recently [17]. 
Interestingly, it has been found that for strongly correlated random variables in different 
types of random walks or fluctuating interfaces extreme value distributions can display sim- 
ple one parameter scaling behavior [lOj [TT] . 

We solve the Master Equation associated with the minimal model and calculate the joint 
probability distribution for C* attaining a maximum value A at time t in the time interval 
[0, T] exactly analytically and semi- analytically. We show that when the system is far from 
the steady state, in the presence of the feedback reaction, the distribution of the maximum 
value A over the time interval [0, T] is spread out over a broader range of A compared to the 
linear model. In contrast, the distribution of the time t when the maximum value occurred 
is much narrowly distributed in the presence of the feedback. This suggests that feedback 
interactions can help single cells to respond sensitively to weak stimulus with a well-defined 
response time even in the face of stochastic fluctuations. 

II. RESULTS 

A. Irreversible kinetics 

In order to understand the role of stochastic fluctuations in affecting the distribution of 
maximum value of C*, it will be instructive to study the deterministic mass action kinetics 
for the concentration of C* (or [C*]) in the reaction, C^C*, which is described by, 

d[C*]/dt = (fci + k p [C*])[C] - k^[CT] • (1) 
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The total concentration, Co = [C] + [C*], is always fixed, and, the rates k\ and k p determine 
timescales for production of (7* from C via a linear first order reaction, and, a second order 
reaction representing a positive feedback, respectively. (7* is converted back to C with a rate, 
k-i. These time scales in a biological network can be regulated by the strength of a stimulus 
that results in generation of (7*, e.g., a weaker (or stronger) stimulus would give rise to longer 
(or shorter) time scales for production of (7*. The rate equation contains a single stable fixed 
point, and, thus starting with any initial concentration, [C*] monotonically reaches a steady 
state determined by the rate constants, and, Co. Consequently, if a reaction initiated with a 
concentration [C*(t = 0)] < [C*(t — >• oo)] is followed until t — T, the maximum value of [C*], 
uniquely determined by the rate constants, Co, and T, is reached at t — T. However, in the 
presence of intrinsic stochastic noise fluctuations, the maximum value of the concentration 
of C* or the time when it is attained will vary in each stochastic 'trajectory' (Fig. 1), where 
every trajectory represents activation of (7* in a single cell. In this situation, P(n,t|m,0), 
the conditional probability of having n number of molecules of the (7* species at any time t 
starting with a distribution P(m, 0) at t = 0, follows the Master Equation, 

dP(n,t\m,0) = ^_ n + ^ + k ^ n _ !))p( n _ 1? t | m? ) + k_ r {n + l)P(n + 1, t\m, 0) 
-{ki{No -n) + k_in + k p n(N - n))P(n, t\m, 0) , (2) 

where, 7V denotes the total number of molecules of C and (7* species. The distribution of 
the maximum number (N) of (7* molecules and the time when the maximum was reached 
in a time interval [0, T] can be calculated by solving of the above Master equation and using 
the renewal equation [20l [21] , 



t 

P(n,t\m,0) = Q N (n,t\m,0) + J dt' F N (t , \m,0)P(n,t\ Nj). (3) 
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Qn{^) t \ m, 0) describes the probability of having n molecules of (7* species at time £, when 
an absorbing boundary condition, Qjsr(n,t\m,0) = for n > TV, is imposed. Fjsr(t\m, 0) 
denotes the probability of arriving at the state n = N for the first time at time t. If the 
time variable is Laplace transformed in the renewal equation, then Fjv(s|m,0) is related 
to P(AT, s|m, 0) simply by, F N (s\m,0) = P(N, s\m,0)/P(N, s\ iV, 0). The unnormalized 
joint probability distribution for attaining a maximum value TV at time t in the time interval 
[0, T] is then given by, 



N 

F N (t\ to, 0)Q N+ i(n, T\ AT, t) . (4) 

n=0 

We then calculate the unnormalized distribution of the maximum value TV over the time 
interval [0, T], i.e., 

Pmax(N,T)= [ dtE N (T,t\m,Q) (5) 
Jo 

and the unnormalized distribution of the time t when the maximum value occurred given 
by, 

Pmax(t,T)= E N (T,t\m,0). (6) 

N=m+1 

We first consider the case where the production of C* occurs irreversibly, i.e., = 0. In 
this limit, E N {T, i, |to, 0) can be evaluated analytically as calculations simplify due to the fol- 
lowing relations: P(N, T\ to, 0) = for m > N, thus, Q N+1 (N, T\ to, 0) = P(N, T\ to, 0) 
for m < N. Consequently, the joint probability distribution can be expressed as, 
E N (T,t,\m,0) = F N (t\m,0)P(N,T\N,t). This essentially implies that the probability 
of having a maximum value TV at time t is the probability the state n = N was attained 
at time t for the first time and then no reaction occurred in the time interval T — t. Next 
we calculate these distributions for the linear and the feedback models by solving Eq. [2] for 
fe_i = 0. 

In the absence of the positive feedback (k p = 0), the exact solution of the 
Master Equation in Eq. § yields, P(N : t\ m, 0)= N °~ m C No - N e ~^ N °~ N ^ klt (l - e - k ^ N - m ^ 
where, P(m, 0) = S m ^. The first passage time distribution is given by, 
F N (t\m) = (N - N + l)k 1 P(N - 1, £| ra, 0), therefore, E N (T,t,\m,Q) = k^N^ - N + 
1)P(N -l,t\m,0)P(N,T\N,t). Thus, P max (N,T) = P(N,T\m,0), and, P max (T,t) = 
(A^o — m)e~ klt (l + e~ klT — e~ klt ) N ° ^ m+1 ^ ( se e [27] for additional details). In the presence of 
the feedback, Eq. [2] can be solved exactly by Laplace transforming the time variable. We 
consider the 'feedback only' {k\ = and k p ^ 0) case to exclusively interrogate the role of 
the positive feedback. The exact solution for the time dependent probability distribution 
for to = 1 is given by, 

N 1 

P (N,4 1,0) = k r(N - i)!"°-'^-„n (8+v(jVo _ r)) m 

The calculation of the inverse Laplace transformation of the above equation is tedious 
but straightforward, and, since the poles of P( iV, s\ 1, 0) at two different values of r 
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can be equal, the probability distribution contains terms which are product of linear 
and exponential functions of t (details in [27J). The first passage time distribution for 
this case is given by, F N (t\m) = k p (N - 1)(N — (N — 1))P(N - l,t\ 1,0), therefore, 
E N (T,t, |1,0) = k p {N - 1)(N -(N- 1))P{N - 1,£| i^Q) e -kpN(N -N)(T-t) m As in the linear 
model we find, P max (A^, T) = P(iV, T\ 1,0). However, P max (t,T) does not possess a simple 
expression as the linear model. The shapes of the distributions, P max (N, T) and P max (£, T), 
depend on N and the dimensionless variable, r = k p T (or k\T for the linear model). In 
order to compare the distributions for the pure feedback and the linear models, we chose an 
end time T, where the average number of C* molecules was the same for both the models. 
Fig. 2a shows the maximum value is distributed more evenly across different numbers of 
molecules of C* in the presence of the feedback compared to the linear model, in contrast, 
Pmax{t->T) (Fig. 2b) is more sharply peaked for the feedback model, indicating that once 
the first molecules of C* are produced the positive feedback leads to fast production of C* 
molecules giving rise to a peak at t — T. The variation (inset, Fig. 2a) of the Fano Factor, 
/ = ((A^ 2 ) — (N) 2 )/ (AT), which quantifies if a distribution is broader than a Poisson distri- 
bution (where, / = 1), with (N) at different times shows that P max (Af, T) is more spread 
out for the pure feedback model as long as the system is away from the steady state. As 
the system approaches the steady state, due to the irreversibility in the reactions, all the C 
molecules are converted into C*, i.e. , P(AT, t — >■ oc|m, 0) — >■ 5tv,tv , and then / decreases and 
become comparable for both the models. We used Kurtosis (K) defined as, K = /a^/a 4 — 3, 
where, /14 and a 2 denote the 4th cumulant and the variance, respectively, to quantify the 
"peakedness" and the presence of "heavy tails" [22] in P max (t,T) compared to a Gaussian 
distribution (K — 0). The pure feedback model produces much larger values of K at different 
values of (N) compared to the linear model (inset, Fig. 2b) indicating higher "peakedness" 
of the distributions in the presence of the feedback. 

The above results become evident at the limit, A^ — >• 00, kiN — >> k\ and k p N — >> k p . 
Then the probability distribution, P(n, t\ m, 0), follows a simple functional form for both 
the models. In the linear model, 

P max (iV, T) = k^-^e-^ihT f^KN - m)\ (8) 

and 

P m ^(T,t) = e-^ T - l \ (9) 
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whereas, in the model with pure feedback[5], 

iW^, T ) = m(m + 1) • • • (N - 1)(1 - e -k P T^N-m e -~k p mT^ N _ m y (10) 

and 

i"max(T,t) = me-^ T e 4 ^ T - f) [l - (1 - e -M) e -^(r-*)]-(m+i) (11) 

Calculation of the Fano factor (/) for P max (N,T) shows that, / = k{T j{m + fciT) < 1 and 
/ = e kpT — 1 for the linear and the feedback model, respectively; clearly demonstrating that 
Pmax(N,T) contains more variation in TV for the feedback model for T > l/k p . On the 
other hand, the denominator in P max (t, T) for the feedback model decreases as t approaches 
T making the distribution more sharpely peaked at t — T compared to the linear model. 
These results provide the following physical understanding. For the feedback model, the 
probability (oc e ~ nkpAt ) to remain in the state n for a time interval At decreases with n, 
whereas, this probability (oc e~ hlAt ) is independent of n for the linear model. Therefore, 
in the feedback model, the states with smaller values of n spend large fraction of their 
time initially in waiting for the first reactions to occur, and then at times closer to the 
end time T, as these states move to larger values of n, the next reactions take place in 
rapid successions. This produces a large range of values of maximum values TV in the time 
interval [0,T] but a sharp peak in P max (t,T). In contrast, in the linear model a state with 
n molecules moves to the next higher state (n + 1 molecules) with a constant rate, so at 
t — T all the states reside close to the average value of (7*. In addition, since the waiting 
time for next reaction to occur does not depend on n, P max (t,T) is more spread out. We 
expect these features to persist even in the presence of a non-zero de-activation rate, as long 
as, the time scales for the feedback reactions are smaller than the de-activation time scale. 
These results suggest the following biological significance. During early signaling events, if 
the rate of de-activation is slower than that of activation, then in time scales smaller than 
that of deactivation the presence positive feedback interactions can help cells to achieve a 
wide range of activation within a narrow response time even in the presence of stochastic 
fluctuations. In the following section we will show the main features of the distributions 
Pmax{N \T) and P max (t,T) persist even in the presence of a non-zero de-activation rate. 
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B. Reversible kinetics 



In the presence of a non vanishing rate of de- activation (i.e., ^ 0), the number of 
C* molecules can decrease after attaining the maximum value, therefore, the simple rela- 
tionships between the P(n,£|m,0) and the maximum value distributions as in the previous 
section no longer hold. Moreover, it becomes difficult to analytically solve the Master Equa- 
tion exactly when the positive feedback is present. Therefore, we calculated the maximum 
value distributions semi- analytically. We briefly outline the method here, the details of the 
calculations are shown in [27J. The Master Equation in Eq.[2] can be cast as an operator 
equation [21] described by, 

d\P(t))/dt = L\P(t)) , (12) 

where, (n\P(t)) = P(n,i|m,0) and (L) nn * = (N — n')(ki + k p n')5 n ^ n _i + Lin'^/^ n+i — 
(ki(N — n) + k-\n + k p nn')5 n ^ n . We solve the above equation by numerically evaluating 
the right (\R r )) and left ((L r \) eigenvectors, and the eigenvalues ({A r }) of the operator, 
L. The solution of the Master equation then can be written as, (n\P(t)) = P(n,£|m,0) = 
(n\e Lt |P(0)) = e Xrt a r (0) (n\R r ) = a r {t)R rn , where, |P(0)) describes the probability 

r=0 r=0 

distribution at t — 0, and, R rn = (n\R r ). {a n (0)} is calculated from the initial condi- 
tion. The same scheme is used to calculate the probability distribution, Qn(^^ t \ m, 0) using 
Qjsr(n,t\m,0) = e x r N)t ar N \o)I$n , where, {A^} and {R^V} are the eigenvalues and 

r=0 

eigenvectors of L with an absorbing boundary condition at n = AT, respectively. We define 

N 

the survival probability, £jv(£| m,0) = Qn(^^ t \ m, 0), which can be used to evaluate the 

71=0 

first passage time distribution, F N (t\m,0) from F/v(i|ra, 0) = —dS^/dt. The maximum 
value distribution functions are then calculated using the following equations, 

T 

Pmax(A^T) = J dtF N {t\m,0)S N+1 (T-t\N,0), (13) 



and, 

Pmax(T,t)= F N {t\m,0)S N+1 (T-t\N,0). (14) 

N=m J +l 

The shape of the probability distributions, P(n, £|m, 0) and E N (t, T|m, 0) depend on the 
dimensionless parameters, k p t ) k\jk v and k_\jk v for the feedback model, and, k\t and k-\jk\ 
for the linear model. We varied k-\ as well as t and T in the models to investigate the effect 
of the de-activation rate on the above distributions. We kept k\jk v fixed to a small non-zero 
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value (0.01) with two goals in mind: (i) prevent the n = state from becoming an absorbing 
state, (ii) to exclusively study the effect of the feedback; when k\/k p 3> 1, the feedback model 
starts behaving like the linear model. The presence of a non-zero de-activation rate makes 
Pmax(N,T) different than P(n, X|ra,0) (Fig. 3a), since P(n,£|m,0) no longer vanishes for 
m > n. The non zero deactivation rate can generate bimodal distributions for P(n, T\m, 0) 
(Fig. 3a) in the feedback model, which is purely generated due to stochastic fluctuations 
since the deterministic rate equation (Eq. [j]) does not possess any bistability. P max (N,T) 
can also show a bimodal distribution as P(n,T|m,0). However, the peaks in P max (N,T) 
occur at larger values of TV compared to the values of n where P(n, T|m, 0) is peaked. This 
effect is produced by stochastic trajectories that attained the maximum value n = N at 
times earlier than T. Since C* in the linear and the feedback models can attain the state 
n = N with a non vanishing probability, P max (iV, T) —> 5n,n as T -> oc in both the 
models, making the distributions similar for both the models at long times. Therefore, 
the bimodal distribution in P max (N,T) for the feedback model will be transient. These 
results can have implications for transient bimodal distributions observed in experiments 
when the underlying signaling network predicts steady state bimodal distributions. In such 
experiments, activation of the cells could be caused by a key signaling protein crossing the 
activation threshold for the first time, consequently, the distribution of activated cells will be 
represented more appropriately by P max (N,T) instead of P(n,T|m,0). Next we calculated 



Pmax{t^T) using Eq. 14 As in the irreversible case, the distribution, P max (t,T), shows 
a higher peakedness for the feedback model as compared to the linear model (Fig. 3b). 
The larger values of the Fano factor (/) for the distribution P max (iV, T) for a range of de- 
activation rates (Fig. 4a and 4c) and end times T demonstrate that P max (N,T) continues 
to have more variance for the feedback model compared to the linear model as long as the 
system is away from the steady state. The peakedness in P max (£,T) was characterized by 
the kurtosis, K ) as in the previous section. The feedback model produces positive and 
substantially larger values of K than that for the linear model (Fig. 4b and 4d), where, for 
most of the parameter values, K is negative, indicating a flatter distribution compared to a 
Gaussian distribution. Therefore, the above results demonstrate that even in the presence 
of de- activation of signaling molecules, the presence of the positive feedback prepares the 
single cells to respond to a wide range of activation thresholds in a well defined time window 
in a noisy environment. 
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III. CONCLUSION 



We have studied how commonly found nonlinear biochemical processes in cell signaling 
and gene regulatory networks such as positive feedbacks influence single cell decision pro- 
cesses in the presence of stochastic fluctuations when these decisions are regulated by key 
proteins attaining threshold concentrations within a time window. We analyzed the joint 
probability distribution of the maximum value of concentration of a molecular species and 
the time when the maximum concentration was reached instead of the probability distri- 
bution of the concentration of the molecular species at any time, as the later distribution 
does not contain information regarding if the molecular species attained the threshold con- 
centration at an earlier time. We calculated the maximum value distributions exactly and 
semi-analytically in minimal models that can effectively describe linear and positive feedback 
interactions in biochemical reactions found in a wide range of cell signaling networks. In 
particular, we investigated the role of positive feedback interactions in affecting the shape of 
the maximum value distributions. We find that in the presence of a positive feedback inter- 
action the maximum values of concentrations of the activated species are distributed more 
broadly compared to the linear model when the system is away from the steady state. How- 
ever, the positive feedback produces a narrower distribution in the time when the maximum 
activation was achieved. Therefore, a positive feedback interaction, even in situations when 
stochastic fluctuations dominate signaling kinetics, provides cells the ability to respond to 
situations when specific cellular proteins need to attain a wide range of threshold concen- 
trations within a narrow time window to influence cell decision processes. This property of 
the positive feedback could play an important role when a cell population has to sensitively 
respond to a weak stimulus. The biological significance of the presence of a "heavy tail" 
at time scales smaller than the most probable time scale in P max (i,T) in the presence of 
positive feedback interactions is less evident. Perhaps, the positive feedback helps create a 
small reservoir of cells that can react to very weak stimuli with a range of relatively smaller 
response times when the majority of the cells that are destined to respond at a much biologi- 
cally irrelevant longer time scale. The probability distributions calculated from the solutions 
of the Master equations for the minimal models indicate presence of multiple time scales 
in the system. It will be interesting to see if this can produce multi scaling behavior [12] 
in extreme value distributions in such chemical reactions in general. Such examples will be 
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qualitatively different than extreme value distributions in standard Brownian motion [H] or 
models of fluctuating interfaces [10] which display single parameter scaling. In addition, the 
minimal models studied here are embedded in larger biological networks which ultimately 
determine cell fate responses, therefore, it will be important to investigate the behavior of 
maximal value distributions when the minimal models are embedded in a larger network. 

This work was funded by the Research Institute at the Nationwide Children's Hospital 
and a grant (1R56AI090115-01A1) from the NIH. I thank C. Jayaprakash and M. Kardar 
for discussions, S. Mukherjee for help with LAPACK routines, and, anonymous reviewers 
for their helpful comments. 
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Figure Captions 



1 1 1 1 1 1 rr 




FIG. 1: Shows three different temporal profiles of the number (n) of C* molecules that attain a maximum value N = 10 in 
the time interval [0, T]. Each profile, representing activation of C* in a single cell, shows cell to cell variation of the time t 
when the maximum value was reached. The stochastic trajectories were obtained by performing a Gillespie simulation 23 26 
of the chemical reaction C C* , where the system started with no C* molecules at t = for the parameters (described in 
Eq.|2l ki = 0.001, k p = 0.01, fc_i = 0.5, and 7V = 50. 
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FIG. 2: Comparison between the feedback and the linear model when k-i = 0, No = 21, and P(m, 0) = <5 m> i. The results 
from the exact solutions (solid lines:feedback model; dashed lines: linear model) are compared with Gillespie simulations 
23 26 of the reactions for the feedback (o) and the linear (□) model. Data from the Gillespie simulations are averaged over 
10 6 stochastic trajectories, (a) Shows P macc (7V, T) at T = 10 for the feedback (k p = 0.01) and the linear (k\ = 0.9699) models 
where both the models produce the same values of (N(T)). (Inset) Variation of the Fano factor / with (N(T)) calculated 
from Pmax(N,T) at different times for the same parameter values as in (a). The feedback model (solid line) displays larger 
values of / as compared to the linear model (dashed line) before the systems reach the steady states, (b) Shows the 
distributions P max (£,T) vs t at T = 10 for the linear and the feedback models. The other parameters are the same as in (a). 
The point at t = shows the probability for the n = 1 state to remain in the same state until t = T. (Inset) Shows variation 
of Kurtosis K with (N(T)) as T is increased. The feedback model displays larger values of K indicating presence of larger 
peaks and heavy tails in Pmax(t,T). 
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FIG. 3: Comparison of the distribution functions P ma x(N,T) and Pmax(T,t) between the feedback and the linear models 
for non zero de-activation rates (k— i 7^ 0). The results are shown for No = 21 and an initial distribution of C* given by 
P(m, 0) = 6m, for both the models, (a) Shows Pmax(N, T) at T = 200 for the feedback (k p = 0.01 , k-i = 0.1 and 
k\ = 0.001) and the linear model 1 = 0.01 and k\ = 0.00269) where both the models produce the same values of (N(T)). 
Data from Gillespie simulations 23 26 averaged over 10 6 stochastic trajectories for the feedback (o) and the linear (□) 
models are compared with the semi-analytical calculations (dark solid lines: feedback; dark dashed lines:linear). P(n, T\m, 0), 
calculated using the semi-analytical scheme for the feedback (gray solid line) and the linear (gray dashed line) models, are 
compared with P m ax(N, T). (b) Variation of P m ax(t,T) with t at T = 200 for the linear and feedback models. The results 
from the semi-analytical calculations are shown in dark solid and dashed lines for the feedback and the linear models, 
respectively. P m ax(0, T) for the feedback model showing the probability of the system to remain at the initial state (m = 0) 
until time t = T is finite and is not shown on the graph. Data from Gillespie simulations, shown using the same visualization 
scheme as in (a), are compared with the semi-analytical calculations. All the parameters used for the calculations are the 
same as in (a). 
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FIG. 4: Comparison of the shapes of the distributions P m ax(N, T) and P max (T,t), characterized by Fano factor, /, and, 
Kurtosis, K, respectively, between the feedback and the linear models for non zero de-activation rates (k—i ^ 0). The results 
shown are calculated using the semi-analytical scheme for No = 21 and an initial distribution of C* given by P(m,0) = 5 mj o. 
The results from the feedback and the linear models are shown with filled circles and filled squares, respectively, (a) Shows 
variation of the Fano factor, /, with (N(T)) calculated using the semi-analytical scheme described in the main text at 
different values of k-i as T is increased uniformly (interval of 50) from 50 to 1000 for the feedback and the linear model. The 
rates used for the feedback model are given by, k p = 0.01, k\ = 0.001, and, k-i = 0.01. The rates for the linear model are set 
at ki = 0.00269 and k-i = 0.0005. The figure shows that the values of / are substantially larger in the feedback model than 
the linear model for the same values of (N(T)) as long as the system is away from the steady state. (Inset) The average value 
of C*, (N(T)), calculated for the same parameters used in (a), increases with T until the system reaches the steady state for 
both the feedback (solid line) and the linear (dashed line) model, (b) Variation of K with (N(T)) for the feedback and the 
linear model are displayed using the same visualization scheme and parameters as in (a). The feedback and the linear models 
produced large +ve, and -ve values of K, respectively. We use log(K) scale for the +ve values of K for better visualization, 
(c) Variation of the Fano factor (/) with (N(T)) for a range of values of k-\. The points depicting variations of / with 
(N(T)) as T is increased at a fixed value of k-i are connected with solid lines. Increasing values of k-i are indicated by 
changing the shade of the symbols from dark to lighter shades of gray, k—i changes uniformly from 0.01 to 0.208 in an interval 
of 0.018 for the feedback. For the linear model, k—i changes evenly from 0.0005 to 0.0995 in an interval of 0.009. All the other 
parameters for both the models are held fixed at the values given in (a), (d) Variation of K with (N(T)) for the feedback and 
the linear model displayed using the same visualization scheme and parameters as in (c). The feedback and the linear models 
produced large +ve, and -ve values of K, respectively. We use log(K) scale for the +ve values of K for better visualization. 



15 



[1] J. Das, M. Ho, J. Zikherman, C. Govern, M. Yang, A. Weiss, A. K. Chakraborty, and J. P. 

Roose, Cell 136, 337 (2009). 
[2] J. E. Ferrell Jr and E. M. Machleder, Science 280, 895 (1998). 
[3] T. S. Gardner, C. R. Cantor, and J. J. Collins, Nature 403, 339 (2000). 
[4] C. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and Natural Sciences, 

(Springer- Verlag, Heidelberg, 2004). 
[5] M. Delbruck, J. Chem. Phys. 8, 120 (1940). 

[6] H. H. McAdams and A. Arkin, Proc. Natl. Acad. Sci. USA 94, 814 (1997). 
[7] L. S. Weinberger and T. Shenk, PloS Biol. 5, e9 (2007). 

[8] P. S. Swain, M. B. Elowitz, and E. D. Siggia, Proc. Natl. Acad. Sci. USA 99, 12795 (2002). 
[9] E. J. Gumbel, Statistics of Extremes, (Dover Publications, Inc, New York, 2004). 
[10] S. N. Majumdar and A. Comtet, Phys. Rev. Lett. 92 225501 (2004). 

[11] S. N. Majumdar, J. Randon-Furling, M. J. Kearney, and M. Yor, J. Phys. A 41 365005 (2008). 
[12] L. P. Kadanoff, Chin. J. Phys. 29 613 (1993). 
[13] J. H. Gillespie, Theor. Popul. Biol. 23, 202 (1983). 

[14] A. Kosmerlj, A. K. Chakraborty, M. Kardar, and E. I. Shakhnovich, Phys. Rev. Lett. 103 
068103 (2009). 

[15] P. Embrechts, C. Kluppelberg, and T. Mikosch, Modelling Extremal Events for Insurance and 

Finance (Springer, Berlin, 2004). 
[16] H. A. Orr, Evolution 56, 1317 (2002). 

[17] M. Artomov, M. Kardar, and A. K. Chakraborty, J. Chem. Phys. 133, 105101 (2010). 
[18] T. W. Burkhard, G. Gyorgyi, N. R. Moloney, and Z. Racz, Phys. Rev. E 76, 041119 (2007). 
[19] S. Redner, A Guide to First-Passage Processes, (Cambridge University Press, Cambridge, 
2001). 

[20] N. Goel and N. Ritcher-Dyn, Stochastic Models in Biology, (Academic Press, New York, 1974). 

[21] J. Honerkamp, Stochastic Dynamical Systems, (VCH Publishers, New York, 1993). 

[22] L. T. DeCarlo, Psychol. Meth. 2, 292 (1997). 

[23] D. T. Gillespie, J. Chem. Phys. 81, 2340 (1977). 

[24] J-P. Bouchaud and M. Mezard, J. Phys. A 30 7997 (1997). 



16 



[25] R. W. Katz and B. G. Brown, Clim. Chan. 21, 289 (1992). 

[26] M. Lis, M. N. Artyomov, S. Devadas, A. K. Chakraborty, Bioinform. 25, 2289(2009). 
[27] See Supplementary Material Document No. for additional details. 



17 



